Outline

• GEOG728 Final Project
• Aerial tourists over National Parks
• Overflights data over Lake Mead, Arizona
  • Flights are polylines
  • Point pattern
• ArcGIS Pro and RStudio

R Packages

library(sf)

## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE

library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.8
## v tidyr   1.2.0     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(pacman)
p_load(tmap)

Lake Mead - Western Part

lake <- st_read("D:/GEOG 728/Project/Lake Western Unit Polygon/LAKE Western Unit Perimeter Polygon.shp")

## Reading layer `LAKE Western Unit Perimeter Polygon' from data source 
##   `D:\GEOG 728\Project\Lake Western Unit Polygon\LAKE Western Unit Perimeter Polygon.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 18 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 686803.9 ymin: 3894772 xmax: 784523.4 ymax: 4052888
## Projected CRS: WGS 84 / UTM zone 11N

glimpse(lake)

## Rows: 2
## Columns: 19
## $ UNIT_CODE  <chr> "LAKE", "LAKE"
## $ GIS_Notes  <chr> "Lands - http://landsnet.nps.gov/tractsnet/documents/LAKE/M~
## $ UNIT_NAME  <chr> "Lake Mead National Recreation Area", "Lake Mead National R~
## $ DATE_EDIT  <date> 2018-05-10, 2018-05-10
## $ STATE      <chr> "NV", "NV"
## $ REGION     <chr> "PW", "PW"
## $ GNIS_ID    <chr> "858265", "858265"
## $ UNIT_TYPE  <chr> "National Recreation Area", "National Recreation Area"
## $ CREATED_BY <chr> "Lands", "Lands"
## $ METADATA   <chr> "https://irma.nps.gov/App/Reference/Profile/2210708", "http~
## $ PARKNAME   <chr> "Lake Mead", "Lake Mead"
## $ Shape_Leng <dbl> 16.6911, 16.6911
## $ GlobalID   <chr> "{375ADCB5-95C2-4E7D-8958-DB04D860493F}", "{375ADCB5-95C2-~
## $ REGS_LINK  <chr> "https://doimspp.sharepoint.com/:w:/r/sites/nps-coast/D_hom~
## $ Shape_Le_1 <dbl> 2078125, 2078125
## $ ORIG_FID   <dbl> 0, 0
## $ Shape_Le_2 <dbl> 994.7563, 806899.9868
## $ Shape_Area <dbl> 6.184385e+04, 6.865387e+09
## $ geometry   <POLYGON [m]> POLYGON ((699441.3 3928336,..., POLYGON ((739767.1 4052888,~

Buffer Lake by 0.5 mile

lake_buff <- st_buffer(lake, dist = 0.5*1609.34)

ggplot()+
  geom_sf(data = lake_buff, aes(fill = UNIT_CODE))

Overflights data - point pattern

flts_pts <- st_read("D:/GEOG 728/Project/LAKE_Mead_Waypoints/LAKE_20200816_to_20210930_Pts_clip_UTM11_comb.shp")

## Reading layer `LAKE_20200816_to_20210930_Pts_clip_UTM11_comb' from data source 
##   `D:\GEOG 728\Project\LAKE_Mead_Waypoints\LAKE_20200816_to_20210930_Pts_clip_UTM11_comb.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 4176388 features and 13 fields
## Geometry type: POINT
## Dimension:     XYZ
## Bounding box:  xmin: 686005.7 ymin: 3893975 xmax: 785269.8 ymax: 4053688
## z_range:       zmin: 2.4e-05 zmax: 6096
## Projected CRS: WGS 84 / UTM zone 11N

dim(flts_pts)

## [1] 4176388      14

glimpse(flts_pts)

## Rows: 4,176,388
## Columns: 14
## $ Operation  <chr> "250820", "250820", "250820", "250820", "250820", "250820",~
## $ Track_Star <chr> "Thu Aug 27 17:18:32 MDT 2020", "Thu Aug 27 17:18:32 MDT 20~
## $ Reg_No_    <chr> "N337NV", "N337NV", "N337NV", "N337NV", "N337NV", "N337NV",~
## $ TYPE_REGIS <dbl> 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,~
## $ NAME       <chr> "SUNRISE ASSET MANAGEMENT LLC", "SUNRISE ASSET MANAGEMENT L~
## $ TYPE_AIRCR <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,~
## $ TYPE_ENGIN <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,~
## $ sinuosity  <dbl> 0.3104282, 0.3104282, 0.3104282, 0.3104282, 0.3104282, 0.31~
## $ Z          <dbl> 878.9993, 880.9706, 886.1369, 881.2022, 873.0740, 871.4991,~
## $ Alt_MSL    <dbl> 2273, 2263, 2256, 2248, 2240, 2233, 2223, 2212, 2199, 2179,~
## $ Alt_AGL    <dbl> 1394, 1382, 1370, 1367, 1367, 1362, 1359, 1347, 1339, 1322,~
## $ Alt_AGL_ft <dbl> 4573, 4534, 4495, 4485, 4485, 4469, 4459, 4419, 4393, 4337,~
## $ Alt_MSL_ft <dbl> 7457, 7425, 7402, 7375, 7349, 7326, 7293, 7257, 7215, 7149,~
## $ geometry   <POINT [m]> POINT Z (690774.6 3986698 2..., POINT Z (690859.9 398~

Point pattern: Altitudes AGL 4000-45000ft

flts_pts_1 <- flts_pts %>%
  filter(Alt_AGL_ft %in% (4000:4500))

dim(flts_pts_1)

## [1] 527584     14

flts_pts_2 <- flts_pts_1[1:10000, ]

dim(flts_pts_2)

## [1] 10000    14

Overflights mapping

ggplot()+
  geom_sf(data = lake_buff, aes(fill = UNIT_CODE))+
  geom_sf(data = flts_pts_2, size=0.005, alpha=0.5)

Mean center

xy <- st_coordinates(flts_pts_2)
mc <- apply(xy, 2, mean)
mc

##         X         Y         Z 
##  711617.0 3990279.6    1794.1

mean_center <- as.data.frame(cbind(mc[1], mc[2]))

ggplot() +
  geom_sf(data = lake_buff) +
  geom_sf(data = flts_pts_2, size=0.05, alpha=0.5) +
  geom_point(data = mean_center, aes(x = V1, y = V2), color = "red", size = 2)

Basemap

tmap_mode("view")

## tmap mode set to interactive viewing

tm_basemap("OpenStreetMap")+
  tm_shape(lake_buff) + tm_polygons("UNIT_CODE", alpha = 0.5)+
  tm_shape(flts_pts_2) + tm_dots(size=0.005, col = "blue")